Report project 05¶

First part¶

  • First we need to build the graphene.
  • Then we need to see how the total energy of the system evolves as we increase the mesh of k-points to sample the reciprocal space.
  • Finally, discuss the parameters that indicate we have achieved convergence.

Second part¶

  • Now we need to build the hexagonal Boron Nitride (hBN).
  • Then plot its bands and compare how they have changed compared to graphene.

To do this exercise you will need:

  • Pseudopotentials for Carbon
  • Pseudopotentials for Boron
  • Pseudopotentials for Nitrogen

The Carbon file will be in the directory Inputs/Graphene

The Boron and Nitrogen files will be in the directory Inputs/hBN

In [1]:
#Import the libraries we will be needing
import numpy as np
import sisl
import sisl.viz
import plotly.express as px
from pathlib import Path
import os
import itertools
import matplotlib.pyplot as plt

Now we will generate the graphene using sisl

In [2]:
# Create the geometry
graphene = sisl.geom.graphene()

Now we will write this structure in a .fdf file

In [3]:
inputs_dir = Path("Inputs/Graphene")
graphene.write(inputs_dir / "graphen_geom.fdf")

Now we define (and create) the directory where we will be studying the effect on changing the k-grid

In [4]:
# Define the root directory for the convergence study for k grid points
root = Path("k_study")
# Create the directory if it doesn't exist
root.mkdir(exist_ok=True)

Now we create a loop that:

  • Creates directories for each size of the k-grid
  • Copies the input files (pseudopotential and structure)
  • Creates a file RUN.fdf where the structure file, the type of basis to generate the wave functions, and the type of k-grid we want to generate in each case are included
  • Collects in a list the total energy of each system
  • Collects in a list the computation time of each simulation
In [7]:
# Set a list with all the k-grid size (in x and y) we want to try
kgrid_points = ['1', '2', '3', '4', '5', '6','7', '8', '9', '12', '15', '18','21', '24']
Es=[]     # Where we will store the total energies of each system
tComp=[]  # Where we will store the Elapsed time

for kgrid_points in kgrid_points:
    # With this part we will create each directory
    k_dir = root / kgrid_points           
    k_dir.mkdir(exist_ok=True)
    
    #print(f"Working on directory with kdrid size of {kgrid_points}") --> for debugging
    
    # With this part we copy the input in each directory
    os.system(f"cp {inputs_dir}/* {k_dir}")
    
    # With this part we create the RUN.fdf file with every part 
    with open(k_dir / 'RUN.fdf', 'w') as f:
        f.write("%include graphen_geom.fdf \n") # Includes the geometry file
        f.write(f"PAO.BasisSize DZP\n")         # Includes the PAO.Basis
        #We now include the kgrid
        f.write(f"%block kgrid.MonkhorstPack \n")
        f.write(f"   {kgrid_points} 0 0  0 \n")
        f.write(f"   0  {kgrid_points}  0 \n")
        f.write(f"   0 0 1  0 \n")
        f.write(f"%endblock kgrid.MonkhorstPack")
    
    # Now we run the simulation
    siesta_path = "/Users/danibedmarromero/anaconda3/envs/siesta/bin/siesta" 
    os.system(f"cd {k_dir}; {siesta_path} RUN.fdf > RUN.out")
    
    # Now we store the energy of the system in the Es list 
    ESystem = sisl.get_sile(k_dir / "RUN.out").read_energy()['total']
    Es.append(ESystem)
    #print(f"Energy of the system= {ESystem} eV")
    
    # Now we store the Elapsed time in the list tComp
    timer =[]
    myfile = open(k_dir / 'RUN.out', 'r')
    myline = myfile.readline()
    while myline:
        if 'Elapsed' in myline:
            timer.append(myline)
        myline = myfile.readline()
    myfile.close()
    
    Timer = timer[0]
    ElTime = float(Timer[39:44:1])
    tComp.append(ElTime)
    
kgrid_points = ['1', '2', '3', '4', '5', '6','7', '8', '9', '12', '15', '18','21','24']
print('Done!')

Gamma-point calculation with interaction between periodic images
Some features might not work optimally:
e.g. DM initialization from atomic data
Job completed
Job completed
Job completed
Job completed
Job completed
Job completed
Job completed
Job completed
Job completed
Job completed
Job completed
Job completed
Job completed

Done!

Job completed

Now we represent how the energy of each system and the computation time of it change according to the size of the k-grid:

In [12]:
figure = px.line(
    x=kgrid_points, y=Es,markers=True
).update_layout(
    title="Convergence of energy and its computational cost with k grid points",
    xaxis_title="k grid points in the x and y directions",
    yaxis_title="Energy [eV]"
)

figure.add_scatter(
    x=kgrid_points, y=tComp, yaxis="y2", line_color="red", name=r"Time [s]"
).add_scatter(
    x=kgrid_points, y=Es, yaxis="y1", line_color="blue", name=r"Energy [eV]"
).update_layout(
    yaxis2 = {"side":'right', 'overlaying': 'y', 'title': "Time [s]"}
).update_layout(
    yaxis = {"side":'left',  'title': "Energy [eV]"}
)

figure.show()

Discussion;¶

We can observe that as the size of the k-grid increases, the total energy evolves towards a constant value (308 eV), which does not change even as the grid size increases. This is a very good indicator that the results have converged because despite increasing our precision, we do not observe considerable changes. In this case, we could say that from a k-grid size of around 7x7 in the x and y directions, we already have enough (even with 4x4 we could be satisfied)!

If we observe the computation time, we notice a very peculiar trend: first, the time decreases until reaching a minimum and then increases again. The second increase in computation time was expected, as an increase in the precision with which we sample the k-space would naturally increase the computational cost. However, the fact that the cost for simulations with low precision is higher than for simulations with an acceptable precision (according to the energy convergence) is surprising. I believe this is due to the system not being able to evolve naturally under the given conditions, leading to an anomalous behavior that is more costly to simulate.

Third part¶

Now we need to simulate the graphene and the hBN to compare them

In [25]:
# Create the hBN by modifying the graphen structure
graphene = sisl.geom.graphene()
hBN = graphene.move([0.03, 0, 0], atoms=1)
  • https://www.researchgate.net/figure/A-BN-hexagonal-primitive-unit-cell-B-BN-Orthorhombic-primitive-unit-cell-C-BN_fig1_309695675
  • https://hal.science/hal-01323899/document
In [26]:
#Write the geometry in a .fdf file
input_dir_hBN = Path("Inputs/hBN")
hBN.write(input_dir_hBN / "hBN_geom.fdf")

input_dir_graphene = Path("Inputs/graphene_study")
graphene.write(input_dir_graphene / "graphene_geom.fdf")

Now we need to modify the file, because we currently have this:

graphene_move.png

And we need to have this:

hBN.png

In [13]:
# Define inputs paths
input_dir_hBN = Path("Inputs/hBN")
# Define the root directory for the hBN study
root_hBN = Path("hBN_study")
# Create the directory if it doesn't exist
root_hBN.mkdir(exist_ok=True)
In [38]:
# Define inputs paths
input_dir_graph = Path("Inputs/graphene_study")
# Define the root directory for the hBN study
root_graph = Path("Graphene_study")
# Create the directory if it doesn't exist
root_graph.mkdir(exist_ok=True)
In [14]:
#Copy the hBN input to where we will run the simulation
os.system(f"cp {input_dir_hBN}/* {root_hBN}")
    
#Creamos el fichero Run.fdf  
with open(root_hBN / 'RUN.fdf', 'w') as f:
    # We include the fdf file that contains the geometry, with the %include statement
    # Note the \n, which means "new line" in text files.
    f.write("%include hBN_geom.fdf \n")
    # We now include the basis size input.
    f.write(f"PAO.BasisSize DZP\n")
    #We now include the kgrid
    f.write(f"%block kgrid.MonkhorstPack \n")
    f.write(f"   7 0 0  0 \n")
    f.write(f"   0 7 0  0 \n")
    f.write(f"   0 0 1  0 \n")
    f.write(f"%endblock kgrid.MonkhorstPack\n")
    f.write("SaveRho true\n")
    f.write("TS.HS.Save true\n")
In [39]:
#Copy the hBN input to where we will run the simulation
os.system(f"cp {input_dir_graph}/* {root_graph}")
    
#Creamos el fichero Run.fdf  
with open(root_graph / 'RUN.fdf', 'w') as f:
    # We include the fdf file that contains the geometry, with the %include statement
    # Note the \n, which means "new line" in text files.
    f.write("%include graphene_geom.fdf \n")
    # We now include the basis size input.
    f.write(f"PAO.BasisSize DZP\n")
    #We now include the kgrid
    f.write(f"%block kgrid.MonkhorstPack \n")
    f.write(f"   7 0 0  0 \n")
    f.write(f"   0 7 0  0 \n")
    f.write(f"   0 0 1  0 \n")
    f.write(f"%endblock kgrid.MonkhorstPack\n")
    f.write("SaveRho true\n")
    f.write("TS.HS.Save true\n")
In [15]:
#Run the simulations
siesta_path = "/Users/danibedmarromero/anaconda3/envs/siesta/bin/siesta" 
os.system(f"cd {root_hBN}; {siesta_path} RUN.fdf > RUN.out")

Job completed
Out[15]:
0
In [40]:
#Run the simulations
siesta_path = "/Users/danibedmarromero/anaconda3/envs/siesta/bin/siesta" 
os.system(f"cd {root_graph}; {siesta_path} RUN.fdf > RUN.out")

Job completed
Out[40]:
0

Proceed with the analysis

In [2]:
RUN_hBN = Path('hBN_study/RUN.FDF')
H_hBN = sisl.get_sile(RUN_hBN).read_hamiltonian()
In [3]:
RUN_graph = Path('Graphene_study/RUN.FDF')
H_graph = sisl.get_sile(RUN_graph).read_hamiltonian()
In [4]:
# We need to define a path of k points
band_struct_hBN = sisl.BandStructure(H_hBN, points=[[0, 0, 0], [2/3, 1/3, 0], [1/2, 0, 0]],
    divisions=100, names=[r"\Gamma", "M", "K"]
)
# Then we can plot the bands
band_struct_hBN.plot()
In [5]:
# We need to define a path of k points
band_struct_graph = sisl.BandStructure(H_graph, points=[[0, 0, 0], [2/3, 1/3, 0], [1/2, 0, 0]],
    divisions=100, names=[r"\Gamma", "M", "K"]
)
# Then we can plot the bands
band_struct_graph.plot()

With the band structure plots, we can appreciate that the main difference (or at least with more direct practical effects on the compound) is that in hBN we find a GAP (of about 4.5 eV), while in graphene we find its famous Dirac cone, which grants it its very special behavior.

In [6]:
# Get the fatbands plot
fatbands_hBN = band_struct_hBN.plot.fatbands()
# Split the contributions by the n and l quantum numbers
fatbands_hBN.split_groups(on="n+l")
In [43]:
# Get the fatbands plot
fatbands_graph = band_struct_graph.plot.fatbands()
# Split the contributions by the n and l quantum numbers
fatbands_graph.split_groups(on="n+l")

With the fatbands, we can observe how the orbitals n=2, l=1 and n=3, l=2, which give rise to the Dirac cone (especially the orbital n=2, l=1), look like. These same orbitals are the ones that contribute to the bands forming the GAP in hBN.

In [19]:
# Get the PDOS plot
pdos_plot_hBN = H_hBN.plot.pdos(
    kgrid=[90,90,1], nE=1000, Erange=[-10, 10],
    distribution={"method": "gaussian", "smearing": 0.1}
)

# Split the contributions by the n and l quantum numbers
pdos_plot_hBN.split_DOS(on="n+l", name="Atom $atoms")
In [44]:
# Get the PDOS plot
pdos_plot_graph = H_graph.plot.pdos(
    kgrid=[90,90,1], nE=1000, Erange=[-10, 10],
    distribution={"method": "gaussian", "smearing": 0.1}
)

# Split the contributions by the n and l quantum numbers
pdos_plot_graph.split_DOS(on="n+l", name="Atom $atoms")

With the pDOS plots, it's very easy to see how in graphene we have orbitals (the p) that contribute to the density of states, whereas in hBN we have a GAP.

In [7]:
rho_hBN = sisl.get_sile(RUN_hBN).read_grid("RHO")
In [8]:
rho_hBN.plot(axes="xy", plot_geom=True)
In [9]:
rho_graph = sisl.get_sile(RUN_graph).read_grid("RHO")
In [10]:
rho_graph.plot(axes="xy", plot_geom=True)

We can appreciate how the electron density is distributed differently, as even though the centers of the hexagons are where the most negative charge density is located, in graphene we have much more charge density on the bonds (thanks to the pi bonds), while in hBN the charge density is completely localized on the Nitrogen atoms. This is because the bonds have a strong donor-acceptor character, as Nitrogen is much more electronegative than Boron.